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Abstract 

We study the continuous absorbing-state phase transition in the one-dimensional diffusive epi- 
demic process via mean-field theory and Monte Carlo simulation. In this model, particles of two 
species (A and B) hop on a lattice and undergo reactions B — > A and A -|- B — > 2B; the total par- 
ticle number is conserved. We formulate the model as a continuous-time Markov process described 
by a master equation. A phase transition between the (absorbing) B-free state and an active state 
is observed as the parameters (reaction and diffusion rates, and total particle density) are varied. 
Mean-field theory reveals a surprising, nonmonotonic dependence of the critical recovery rate on 
the diffusion rate of B particles. A computational realization of the process that is faithful to the 
transition rates defining the model is devised, allowing for direct comparison with theory. Using 
the quasi-stationary simulation method we determine the order parameter and the survival time 
in systems of up to 4000 sites. Due to strong finite-size effects, the results converge only for large 
system sizes. We find no evidence for a discontinuous transition. Our results are consistent with 
the existence of three distinct universality classes, depending on whether A particles diffusive more 
rapidly, less rapidly, or at the same rate as B particles. 
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I. INTRODUCTION 

This work is devoted to the one-dimensional diffusive epidemic process (DEP) a 
model system in which two kinds of particles, A and B, diffuse on a lattice and undergo 
reactions B — A and A + B — > 2B. There is no intrinsic limit on the number of particles that 
may be present at a given site; the total number of particles is conserved. In the epidemic 
interpretation A represents a healthy organism and B and infected one, with the reactions 
above corresponding, respectively, to spontaneous recovery and transmission of disease on 
contact. Other interpretations are possible, for example, A could represent a properly folded 
protein and B a misfolded one, etc. 

The DEP is a nonequilibrium stochastic model exhibiting a phase transition to an absorb- 
ing state 0,0,0,0]. Such phase transitions arise in many models of epidemics, population 
dynamics and autocatalytic chemical reactions, and have attracted much interest in nonequi- 
librium statistical mechanics, in efforts to characterize the associated universality classes. 
The simplest example is the contact process (CP), or its discrete-time version, directed per- 
colation (DP). In the CP, each site of a lattice is either vacant (0) or occupied by a particle 
(X). Particles die (i.e., the reaction X — 0) at rate 1, independent of the configuration of 
the rest of the system, and reproduce at rate A (reaction X + — 2X) . The offspring parti- 
cle survives if and only if the site it is sent to (selected at random from the neighbors of the 
reproducing particle) is vacant. The CP is therefore a minimal model for birth-and-death 
processes with local competition for space. The particle-free configuration is absorbing. It 
is known that for reproduction rate A < Ac, the stationary density of particles p is zero, 
and that (in the infinite system size limit), p grows continuously from zero as A is increased 
beyond Ac. Thus p serves as the order parameter for this phase transition. 

Critical scaling in the contact process and allied models has been studied extensively, 
both theoretically and numerically. A central conclusion deriving from these studies is 
that critical behavior of the DP type is generic for models exhibiting a continuous phase 
transition to an absorbing state, in the absence of any additional symmetries or conserved 
quantities. (Note that the absorbing state may in some cases encompass more that one 
configuration, as is the case in the pair contact process jo^. If these configurations are not 
related by any symmetry the critical behavior is still expected to fall in the DP universality 
class.) Particularly relevant in the present context is the fact that the diffusive CP (in which 



particles hop at finite rate) also belongs to the DP class 0]. 

The DEP offers a more complicated scenario of scaling. Let Da {Db) denote the diffusion 
rate of A (B) particles. An especially interesting aspect of the DEP is that the absorbing- 
state phase transition appears to belong to three distinct classes, depending on whether 
Da < Db, Da = Db, or Da > Db- The renormalization group (RG) analysis of Oerding et 
al. ^ predicts a continuous phase transition in the first two cases. But for Da > Db there 
is no fixed point, leading to the conjecture that the transition is discontinuous in this case. 
These authors provide numerical evidence of such a transition in two dimensions; numerical 
studies of the one-dimensional DEP, however, show a continuous phase transition l9|. 

The diffusive epidemic process was initially studied via RG methods . Numerical 

simulations for equal diffusion rates were reported in Ref. |u|, yielding results in disagree- 
ment with the RG prediction = 2 (see Subsequently, the simulations reported 
in appeared to resolve this point, but suggested other departures from RG predictions in 
the one-dimensional case. 

Given the disagreement between theory and simulation, it is of interest to perform further 
analysis of the DEP. In the present work we study a one-dimensional version of the process, 
formulated in continuous time, in a manner allowing for a simple formulation of the master 
equation. We study this process (defined in detail in the following section), using one- 
and two-site mean-field theory (Sec. Ill), and Monte Carlo simulation (Sec. IV) using the 
quasi- stationary simulation method. Conclusions and open questions are discussed in Sec. 
V. 



II. MODEL 

The DEP is defined on a lattice of L'^ sites. A configuration is specified by the set of 
variables aj and bj, denoting the number of A and B particles at each site j. There are 
no intrinsic restrictions on the number of particles at each site, making this process, in the 
language commonly employed in the literature, a "bosonic" model. (To avoid confusion we 
note that quantum statistics do not enter into the problem.) The model is a continuous-time 
Markov process characterized by four kinds of transitions: 

• Hopping of A particles to a randomly chosen nearest neighbor (NN) site, at rate Da- 



• Hopping of B particles to a randomly chosen NN site, at rate Db- 

• Transformation of B particles to A particles, at rate r. 

• Transformation of A particles to B particles, in the presence of a B particle at the 
same site, at a rate of A per A-B pair. 

This means that a given site j loses (via diffusion) an A particle at rate Dj^aj (similarly for 
loss of a B particle), undergoes the process B ^ A at rate rbj, and the process A + B — >• 2B 
at rate ^cijhj. Note that all transitions conserve the total particle number = + 

The process involves a rather large set of parameters: Da, Db-, r, A, and the particle 
density p = N/L'^. Since one of the rates may be eliminated through a suitable scaling of 
time, we set A = 1 from here on. This still leaves four control parameters. In the studies 
described below, we fix the diffusion constants and study the behavior in the r-p plane, or 
fix Da and p and treat Db and r as the control parameters. 

III. MEAN-FIELD THEORY 

We begin by studying the model using dynamic mean-field theory (MFT) P]. Here we 
present the results for one- and two-site approximations in one dimension, which follow 
from the master equation for the probability distribution P{ai,bi] aL,bL,t). Consider the 
evolution of the one-site marginal distribution P{a,b): 

P{a, b) = D^iia + l)P{a + 1, 6) - aP{a, b)] + DaY^ a'[P{a - 1, b; a', b') - P{a, b; a', b')] 

a',b' 

+ DB[{b + l)P(a, 6 + 1) - bP{a, b)] + DbJ2 b'[P{a, b - 1; a', b') - P{a, b; a', b')] 

a',b' 

+ r[{b+l)P{a-l,b+l)-bP{a,b)] + {b~l){a + l)P{a+l,b-l)-abP{a,b) (1) 

Here P{a, b; a', b') is the joint probability distribution for a pair of NN sites. Eq. is the 
first in a hierarchy of equations for the probability distributions of 1, 2,...,n,... sites. At 
each level of the hierarchy, the diffusion terms couple the n-site distribution to that for 
n + 1 sites. One-site MFT truncates this hierarchy at lowest order, via the factorization 
P(a, 6; a', b') = P{a, b)P{a', b'), leading to. 



P(a, b) = DA[{a + l)P{a + 1,6)- aP(a, b)] + DApA[P{a - 1, 6) - P{a, b)] 
+ DB[{b + l)P(a, 6 + 1) - 6P(a, b)] + DBpB[P(a, 6 - 1) - P(a, b)] 
+ r[{b+l)P{a-l,b+l)-bP{a,b)] + {a + l){b-l)P{a + l,b-l)-abP{a,b) (2) 

where pA = J2a b (^P{(^, b) is the density of A particles and similarly for pB- An equation for 
Pa is found by multiplying the above equation by a and summing over all values of a and b, 
giving, 

Pa = rps - (ab) (3) 

where (a^fo") = Xla 6 '^'"^"-^('^' (Thus (a) = and similarly for pB-) Note that, as 
discussed in detail below, the cross-moment {ab) is in general different from the simple 
product of A and B particle densities. If we nevertheless set {ab) = paPb, we have 

Pa = rpB - pApB (4) 
With the constraint Pa + Pb = Pi constant, we then find 

pB = {p- r)pB - Pb (5) 

showing that at this (lowest) level of approximation the order parameter pB satisfies the 
Malthus-Verhulst equation with reproduction rate p — r. At this level, which may be called 
a "rate equation" , a continuous phase transition occurs at p = r, independent of the diffusion 
rates; the stationary density of B particles follows 'Pb = P — r. While the more detailed 
mean-field approximations described below yield a more reliable prediction for the phase 
boundary, we observe, in all continuous phase transition, and a linear relation 

between and p — Tc in the vicinity of the critical point (i.e., the usual mean- field critical 
exponent j3 = 1). We find no hint of a discontinuous transition. 

A somewhat better result is obtained integrating (numerically) the full set of one-site MF 
equations, Eq. Q. In numerical analysis, it is necessary to truncate this set of equations at 
cutoff values Oc and be- This is justified since the probability distribution falls off exponen- 
tially for large a and/or b. The cutoff leads to certain technical restrictions in the numerical 
analysis. Naturally, transitions of the form a^ —>■ ac + 1 must be excluded from consideration. 
Moreover, a transition of the form a — a — 1, due to an A particle hopping away from the 



site of interest, must have its rate multiplied by 1 — P4(ac), where PaI*^) = J2b -^('^' ^) 
one-site marginal distribution for the number of A particles. (Similar restrictions apply to 
transitions involving B particles.) 

Once we take the full one-site probability distribution into account, the results for the 
phase boundary Pc{r) depend on the diffusion rates. Two factors enter into this dependence. 
First, in the vicinity of the phase transition, the reaction terms cause the marginal distribu- 
tion for the number of B particles to deviate significantly from a Poisson distribution, while 
rapid diffusion tends to make the distribution more Poisson-like. Second, the reactions cause 
the variables a and b to be anti-correlated [that is, cov(a, 6) = (ab) — PaPb < 0], whereas 
rapid diffusion tends to eliminate this correlation. 

Figure Q] shows the critical line Pc{r) as predicted by the one-site MFT for various com- 
binations of Da and Db- The higher the diffusion rates, the more closely pc approaches the 
simple rate equation result pc = r. For finite diffusion rates pc is always greater than p, due 
again to the anti-correlation of a and b. The critical value pc appears to be more sensitive 
to Db than Da- In fact, for Db = 0.2, the curves for Da = 0.2 and Da = 1 are virtually 
identical. (For larger values of Db, increasing Da does reduce the critical density.) 




MFT1 

FIG. 1: Mean-field predictions for the critical particle density pc versus recovery rate re- Solid lines, top to 
bottom: Da — Db ~ 0.2; Da ~ 0.2, Db = 1; Da = 1, Db = 1; Da = Db = 5; rate equation result, pc = r. 
Bold dashed line: two-site mean-field theory for the case Da — I, Db = I- 



A richer description is obtained when we extend the approximation to two sites. 
There arc (in general) 16 transitions into (and out of) a given state {a,b; a' ,b'). Using 
the symmetry P{a,b;a',b') = P{a',b';a,b), and factoring three-site probabihties so that 
P(a, b; a', b'; a", b") — P{a, b; a', b')P{a', b'; a", b")/P{a', b'), the equation governing the two- 
site joint probabihty can be written as 

dP{aM a',b') ^ ^ ^ ^ ^ ^. y ^ , ^ ^. ^, ^ y 

at 1 

+(a' + l)P(a - 1, b] a' + 1, b') + (a + l)P(a + 1, 6; a' - 1, b')] 

^[{b + l)P(a, b + 1; a', b') + (6' + l)P(a, 6; a', 6' + 1) 

+ {b' + l)P(a, 6 - 1; a', 6' + 1) + (6 + l)P(a, b + 1; a', 6' - 1)] 

+^ [$A(a - 1, &)^(a - 1, 6; a', b') + $^(0' - 1, a' - 1, &')] 

+^ [$B(a, 6 - l)P(a, b - 1; a', 6') + ^^(a', 6' - l)P(a, b- a', 6' - 1)] 

+r[(6 + l)P(a -1,6+1; a', 6') + (6' + l)P(a, 6; a' -1,6' + 1)] 

+ (a+l)(6-l)P(a + l,6-l;a', b') + (a'+ 1)(6'- l)P(a, 6, a'+ 1, 6'- 1) 

- ^DA{a + a') + DB{b + 6') + ^[^^(a, b) + $^(a', 6')] 

+ ^[$B(a, 6) + $s(a', 6')] + r{b + 6') + a6 + a'6'| P(a, 6; a', 6') (6) 

where 

is the conditional A-particle density at a site, given that one of its nearest neighbors has 
occupancy (a, 6). ($5 is defined analogously.) The above equations are integrated numer- 
ically using the fourth-order Runge-Kutta method and a cutoff of 10 for the variables a, 
b, a' and b'. For densities p < 2, the error incurred is negligible. Figure 1 shows that the 
two-site approximation predicts a larger value of pc than does the one-site approximation, 
other parameters being equal. 

We compare the mean-field predictions for the critical recovery rate against simulation 
for three representative cases (all for density p = 1) in Table I. The site approximation 
overestimates Vc by a factor of up to 3.3; in each case the two-site approximation yields a 
substantial improvement, although it still overestimates the critical value by a factor of up to 
2.3. In principle, further improvement could be furnished using higher order approximations. 



but in the present case the computational demands seem excessive. [The number of equations 
to be integrated in the ra-site approximation is [(oc + l)(&c + I)]"-] 

An important aspect of the model is the anti-correlation between variables a and 6 at a 
given site. To quantify this we study 



At the critical point {ps —>■ 0), Qa represents the excess density of A particles at a site 
bearing a B particle. In the one-site approximation, we find (at the respective critical 
points), Qa = -0.583, -0.298, and -0.0895, for Da = Db = 0.2, 1, and 5, respectively. 
That is, the species are anti-correlated, and the magnitude of the correlation decreases with 
increasing diffusion rate, as expected. This effect appears even more markedly in the two-site 
approximation, where, for example, Qa = —0.447 at the critical point with Da = Db = ^■ 
(All results quoted are for density p = 1.) Similar values are found in simulations. The 
two-site approximation shows that the anti-correlation of the numbers of A and B particles 
extends to the nearest neighbor site: {ajbj^i) / pB — Pa = —0.444 for the parameters noted 
above. The variables bj and on the other hand, show a strong positive correlation. 

Since diffusion of B particles is essential to the survival of the process, one might suppose 
that as is reduced, the critical density would increase (for fixed recovery rate r), or 
the critical recovery rate decrease (for fixed density). Mean- field theory (at both levels) 
predicts otherwise, as shown in Fig. |2l the critical recovery rate exhibits a minimum at a 
small value, D^, of the B diffusion rate, but for even smaller values, it grows rapidly. For 
Db > D* the critical recovery rate grows systematically, and appears to saturate at the 
rate-equation value, = p. Fig. El shows that follows the same qualitative trends in 
simulations as in the mean-field approximations, although its numerical value is (generally) 
considerably smaller. 

The reason for the increase in rc at small values of Db appears to lie in the possibility 
of accumulating many B particles at a single site. Suppose that, by some fluctuation, a site 
acquires several B particles. Since Db Da, this accumulation is relatively long-lived, and 
any A particles at this site will rapidly change to B, since the effective rate for the reaction 
A ^ B, at this site, is 6, which is large compared to r, the rate of the inverse reaction, 
and also large compared to Da- (Note as well that A particles straying onto such a site will 
readily change to B, and tend to remain at this site for a long time.) Indeed, for Db < D*, 



Qa = 




- PA 




(8) 
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FIG. 2: Critical recovery rate versus B particle diffusion rate Z?b, for Da — 0.5 and p — 1- Upper curve: 
one-site MFT; lower curve: two-site MFT; points: simulation. 

mean-field theory reveals that the one-site marginal distribution P(&), though quite small for 
6 > 1 decays very slowly with increasing h. The distributions P{a) and P{h) are compared 
in Fig. El (for the case p = 1, r ~ Tc = 0.688, Da = 0.5 and Db = 0.02). Although (globally) 
almost all particles belong to species A, P{h) decays much more slowly than P{a) for large 
occupancies. (In this regime we increase Oc and he to 40 in the one-site MF calculations, to 
ensure that the cutoff does not affect the result.) P{a) is well approximated by a Poisson 
distribution, while P{h) shows strong deviations from this from. As shown in Fig. El (inset), 
similar distributions are observed in simulation, in the regime -D_b -C Da- 



TV. SIMULATIONS 



We perform Monte Carlo simulations using a simulation algorithm designed to reproduce 
faithfully the transition rates defining the process. Simpler, more efficient computational 
models involving the four types of transition [A and B diffusion, recovery (R), and infection 
(I)], are possible, but do not correspond to the same set of transition rates. Our simulation 
method permits quantitative comparison with theoretical predictions, including systematic 
expansions of the master equation jl4l] . 




FIG. 3: One-site probability distributions for the number of A and B particles. Main graph: one-site MFT, 
p = 1, r ~ rc — 0.688, Da = 0.5 and Db = 0.02. Inset: simulation result for r = 0.5 and the same values of 
p, Da and Z?b, system size L = 500. 

The simulation consists of a sequence of events, each of which involves choosing the type 
of transition and the site at which it occurs. The choice of event type depends the total 
transition rates for each of the four processes, given by: 

• Hopping of A particles: total transition rate Wa = NaDa, where Na = % is the 
total number of A particles. 

• Hopping of B particles: total transition rate Wb = NbDb- 

• Transformation of B particles to A particles: total transition rate = tNb- 

• Transformation of A particles to B particles: total transition rate Wi = Y2j (^j^j- 

If we let Wt = Wa + Wb + Wr + Wi denote the total transition rate for all processes, 
then the probability that the next event is of type m (= A, B, R or I), is Pm = Wm/Wr, 
while the mean time to the next event is At = 1/Wt- (After each event the time is advanced 
by this amount.) The next event is chosen at random from this set of probabilities. Once 
the event type is determined, the site at which it occurs must be chosen. For this purpose a 
number of lists are maintained. For example, if the chosen event is hopping of an A particle. 



we select a site at random from a list of all sites (call it the A-list) having aj > 0. Not all 
sites on the A-list, however, are equally likely to host the event: since each A particle has 
the same hopping rate, the probability of the event occurring at site j is proportional to 
ttj. We therefore keep track of the number a^ax of A particles at the site with the largest 
number of such particles. When site j is chosen from the A-list, it is accepted for the next 
event with probability p^cc = o-j/o-max- In case of rejection, we again select a site (/c, say) 
from the A-list and compare with a^ax- In his manner we ensure that each A particle 
has the same likelihood of hopping. In a diffusion event the target site is chosen at random 
from the nearest neighbors of the original site. 

The same procedure is adopted for the other three processes, necessitating the mainte- 
nance of a B-list and an AB-list (the latter containing all sites such that ajhj > 0). The 
rejection procedure outlined above, while necessary to maintain fidelity to the original tran- 
sition rates, is computationally expensive. For this reason we restrict the present study to a 
relatively low density, p = 1, since in this case the maximum values, amax, bmax, and {ah)max 
are generally not very large. (One could in fact use lists of the positions of all A and B 
particles, instead of A- and B-bearing sites. But since the A-B reaction step requires an 
AB site list, we adopt a uniform procedure for all processes. Using particle instead of site 
lists for all processes except the A-B reaction could improve efficiency, particularly for large 
diffusion rates.) 

In the studies reported here we sample the quasi- stationary (QS) distribution of the 
process, (that is, conditioned on survival), an approach that has been found very useful 
in the study of systems with an absorbing state [l5|. (In fact, conventional simulations of 
"stationary" properties of models with an absorbing state actually study the quasi-stationary 
regime, given that the only true stationary state for a finite system is the absorbing one). We 
emp loy a recently devised simulation method that yields quasi-stationary properties directly 
jig . This is done by maintaining, and gradually updating, a set of configurations visited 
during the evolution; when a transition to the absorbing state is imminent the system is 
instead placed in one of the saved configurations. Otherwise the evolution is identical to 
that of a conventional simulation. (The set of saved configurations is updated by replacing 
one of the saved configurations with the current one, with a small probability prep at each 
time step.) 

n 

The above scheme was shown [16|| to yield precise results, in accord with the exact QS 



distribution for the contact process on a complete graph, and with conventional simulations 
of the same model on a ring jisl. The scheme has also been shown to yield results that 
agree, to within uncertainty, with the corresponding results of conventional simulations for 



a sandpile model 17[ . The advantage of the method is that a realization of the process can 
be run to arbitrarily long times. Thus, whereas in conventional simulations a large number 
of realizations must be performed to have a decent sampling of the quasi-stationary state, 
here every realization provides useful information. This leads to an order of magnitude 



3- 



improvement in efficiency in the critical region. For further details on the method see 

We performed extensive simulations of the DEP on rings of L = 200, 500, 1000, 2000, and 
4000 sites, using the QS method. The number of saved configurations M ranges from 1000, 
for L = 200, to 100 for L = 4000. Values of Prep range from 10 ^ to 10 ^ (smaller values 
for larger systems). Two time scales appear to be relevant in the choice of the replacement 
probability Prep- One is the mean residence time of a configuration on the list, = M/p^ep- 
The other is the QS lifetime r, i.e., the mean time between attempts to visit the absorbing 
state. We find that our results are independent of the choice of Prep provided r/r^ < 1. This 
appears to be associated with the need to preserve configurations visited prior to the last 
attempted transition to the absorbing state. Of course one could make tl arbitrarily large 
by reducing Prep, but this would prolong the memory of the initial state. (For this reason, 
we use a Prep ten times larger during the relaxation phase. The latter represents about 10% 
of the total simulation time.) 

Initially, half the particles are of type A and half B; the particles are distributed randomly 
and independently over the sites, so that the distribution of a and 6 at a given site is essen- 
tially Poisson with mean 1/2. Each realization of the process runs for a certain maximum 
time, Tm, of up to 10® time units. The results reported here represent averages over 4 - 
12 independent realizations for each set of parameter values. Averages are taken in the QS 
regime, after discarding an initial transient, with a duration that depends on the system 
size. In practice we accumulate histograms of the time during which the system has exactly 



1, 2,...n,..., B particles. T 
Pb and the moment ratio 



le histograms are used to evaluate the mean B particle density. 
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m=^^ (9) 
(6)2- ^""^ 

We determine the critical recovery rate r^ip, Da, Db) using the criteria of power-law 



dependence of and r on system size L (i.e., the usual finite-size scaling relations ~ 
]^-l3/i^± and T L^). To probe the three characteristic regimes of the process, three cases 
are studied in detail: Da = Db = 0.5; Da = 0.5, Db = 0.25; and Da = 0.25, Db = 0.5. 

Consider first the case Da = 0.5, Db = 0.25. Since it has been suggested that the 
phase transition is discontinuous for Da > Db, it is of interest to study the behavior of 
the order parameter pb as a function of the recovery rate r. In Fig. |3]we plot estimates 
for the order parameter in the limit L oo, based on our data for L = 200 - 4000. (The 
extrapolation can only be performed reliably for r < 0.22 with the data at hand.) It appears 
that the order parameter decays continuously to zero as r ^ ~ 0.23. (We cannot rule 
out a weakly discontinuous transition on the basis of these data, but a continuous transition 
seems the more natural interpretation.) The evidence for a continuous transition is greatly 
strengthened by our observation of critical scahng, as we now discuss. 
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FIG. 4: Order parameter ps versus recovery rate for Da = 0.5, Db — 0.25, and p = 1- The points represent 
extrapolations to the infinite-size hmit based on data for sizes L — 200 - 4000. 

Data for the QS order parameter pB versus system size L are shown on log scales (for 
the case Da = 0.5, Db = 0.25) in Fig. El The data for r = 0.231 curve upward, while those 
for r = 0.234 curve downward, leading to the estimate Tc = 0.2325(10). For this range of 
values we see good evidence of power-law scaling, as expected at a critical point. Similarly, 
the data for the QS lifetime r (see Fig. ^ show power-law scaling for r ~ r^, and significant 



curvature away from Vc- 




FIG. 5: Order parameter pB versus system size L, for the same parameters as Fig. 4. From top to bottom: 
recovery rate r = 0.231, 0.232, 0.233, 0.234. 




FIG. 6: Lifetime r versus system size L, symbols and parameters as in Fig. 5. 

As a check on our procedure for determining rc, we also perform, for Da = 0.5, Db = 0.25, 
initial decay studies [3|2^. In these studies the order parameter ps is followed as a function 



of time, using the same initial condition as in the QS studies. (Due to the large system size, 
the process does not enter the absorbing state on the time scale of the simulation.) Here 
the expected critical behavior is ~ t~^. Using deviations from the power law to identify 
off-critical values, a study of systems of 10^ sites, to a maximum time of 10^, enables us to 
restrict Vc to the interval [0.231, 0.239]. This is consistent with, but less precise than, the 
QS results. Analysis of the data for t > 5 x 10^ furnishes 6 ^ 0.5. The lack of precision in 
these studies and in our result for the moment ratio highlights the slow convergence of these 
simulations, presumably reflecting strong finite-size effects. 

The studies described above lead to estimates for the critical exponent ratios /3/z/_l and 
z = I'w/v.L- To obtain a direct estimate of z/^, we study the order parameter in the neigh- 
borhood of the critical point for various lattice sizes. This permits us to estimate the 
correlation length exponent z/^, using the finite-size scaling relation for the order parameter, 



where A = r — Tc and JF is a scaling function. This implies 



(9 In pb 



dr 



oc 



(10) 



(11) 



For the case Da = 0.5, Db = 0.25, the data for L = 500, 1000, 2000 and 4000 yield the 
estimate z/_l = 2.3(3). We performed similar finite-size scaling analyses for the other cases. 

The moment ratio m is also useful in characterizing critical behavior. Our results, shown 
in Fig. 13 lead to estimates for the hmiting (L — > oo) value except in the case Da = 0.5, 
Db = 0.25, for which the moment ratio decreases systematically with system size L. The 
trend to smaller values of m with increasing Db/Da may reflect the reduced tendency 
for B particles to accumulate at a given site. Our results for the critical parameters are 
summarized in Table II. 

We close this section with the observation that the DEP is characterized by unusually 
large and long-lived fluctuations, which make it difficult to obtain precise results in simula- 
tions. In the case of the contact process, for example, a system size of L = 1000 is sufficient 
for the determination of critical exponents to rather high precision. As in other models dom- 
inated by diffusion, such as conserved sandpiles j^| and the diffusive pair contact process 
(PCPD) 0, 22], scaling properties of the DEP emerge clearly only at larger sizes. One 
signal of this is the slow convergence (with L) of the order parameter moment ratio m. 
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FIG. 7: Moment ratio m versus system size L, at the critical point rc, for (upper to lower) Da — 0.5, 
Db = 0.25; Da^Db^ 0.5; and Da = 0.25, Db = 0.5. 

V. DISCUSSION 



We study the one- dimensional diffusive epidemic process using mean-field approximations 
and Monte Carlo simulations. Mean-field theory provides a qualitative description of the 
phase diagram, and reveals a surprising, non-monotonic dependence of the critical recovery 
rate on the B particle diffusion rate, confirmed in simulation. This reflects the tendency 
(for Db Da) for many B particles to accumulate at a single site, as revealed in the 
one-site probability distribution P{h). Mean-field theory also captures qualitatively the 
anti-correlation in the number of A and B particles at a single site. 

Although our simulations do not extend to sufficiently large systems to furnish precise 
values of critical exponents, our results, especially for the ratio clearly support the 

scenario of distinct universality classes for the three cases, D^ < Db-, Da = Db-, and 
Da > Db- The transition appears to be continuous in all cases. Renormalization group 
analysis [l|, Is, 10, 1^ predicts (independent of the relative magnitudes of Da and Db) that 



z = 2 and 1/^ = 2/(1. The simulation results are consistent with these predictions except ior 
the case Da = 0.25, Db = 0.5, for which our estimates are significantly smaller. This may 
be due to the slow convergence noted above. We hope to address this point in future studies 
of larger systems. It would also be of interest to verify that the critical exponents are in fact 



insensitive to changes in the diffusion rates, within the three universality regimes that have 
been established. We expect the present results to serve as a point of reference for studies 
based on systematic expansions of the master equation. 
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Table I. Critical recovery rate Tc in one- and two-site 

mean-field approximations, compared with simulation. 



Da Db Tc (1-site) Vc (2-site) Vc (sim) 



0.5 0.25 0.5420 
0.5 0.5 0.5771 
0.25 0.5 0.5144 



0.411 0.2325(10) 
0.429 0.1921(5) 
0.368 0.1585(3) 



Table II. Critical parameters from simulation. 

D y Db - in 

0.5 0.25 0.404(10) 2.01(4) 2.3(3) < 1.15 

0.5 0.5 0.192(4) 2.02(4) 2.0(2) 1.093(10) 

0.25 0.5 0.113(8) 1.6(2) 1.77(3) 1.06(1) 



